##################################################################################################################################################################################################################
# Custom Functions
##################################################################################################################################################################################################################

# Block Bootstrap  
blockboot <- function(df, sims = 1000, model, clustername, status, wt = "noweight", inter = "none"){
  holdbeta <- NULL
  
  for(s in 1:sims){
    # Choose clusters to sample
    
    # Make sure we get correct number of preclearance and non-preclearance states -
    # there are going to be some states that have preclearance counties within them but these might be sampled in the uncovered states, with the covered counties still marked as 1
    
    comp <- data.frame(table(df[, clustername], status))
    
    # Identify States in Treatment and Control Groups; Ultimately We Want Balance Here  
    t.clusters <- as.character(comp$Var1[comp$status == 1 & comp$Freq !=0])
    c.clusters <- as.character(unique(comp$Var1[!(comp$Var1 %in% t.clusters)]))
    
    # Get a Sample with Repeats  
    samp.clusters <- data.frame(state_name = as.character(sort(c(sample(t.clusters, length(t.clusters), replace = T), sample(c.clusters, length(c.clusters), replace = T)))))
    
    # Grab this sample out of the original data with repeats  
    bootdat <- apply(samp.clusters, 1, function(x) {df[which(df[ , clustername] == x),]})
    
    # Give every state a unique state ID. The reason to do this is that we want unique clusters. If I draw 3 Texases, for instance, I want 3 separate clusters since these are meant to be 
    # 3 distinct states despite having the same data
    for(i in 1:length(bootdat)){
      bootdat[[i]]$unitid <- paste(i, bootdat[[i]]$fips, sep = ".")
    }
    # Note that this generates a "unitid" field in the function. The model form is going to call a unitid variable from within the bootdat, so you'll want to have something
    #by this name in the bootstrapped data. This, as well as the samp.clusters line, also requires your df to have a column called "State_Name". We can always make the function more general to
    # let you specify what the cluster is, but I'll add that later if we need it.
    
    # Put the list into a single data frame  
    bootdat <- bootdat %>% bind_rows()
    
    # Run Regression on Data Using FELM  
    mod.form <- sub(".*= *(.*?) *=.*", "\\1", model$call)[2]
    if(wt == "noweight"){
      bootholder <- felm(as.formula(mod.form), data = bootdat)
    } else {
      bootholder <- felm(as.formula(mod.form), weights = bootdat[, wt], data = bootdat)
    }
    
    holdbeta <- rbind(holdbeta, bootholder$coefficients)
    print(s)
  }
  # Report Results
  # These are hardcoded to return coefficients for the interactions on preclearance and shelby. If your interaction is called something else, this will fail.
  boot.res <- data.frame(vars = colnames(holdbeta), ses = apply(as.matrix(holdbeta[row.names(holdbeta) == inter, ]), 2, sd, na.rm = T), row.names = NULL)
  if(inter == "none"){
    boot.res$orig.coef <- coef(model)["preclearance:vra"]
    boot.res$orig.se <- summary(model)$coefficients["preclearance:vra" , 2]
  } else {
    boot.res$orig.coef <- coef(model)[inter]
    boot.res$orig.se <- summary(model)$coefficients[inter , 2]
  }
  boot.res$lbounds <- boot.res$orig.coef - qnorm(0.975) * boot.res$ses
  boot.res$ubounds <- boot.res$orig.coef + qnorm(0.975) * boot.res$ses
  boot.res$zs <- abs(boot.res$orig.coef / boot.res$ses)
  boot.res$ps <- 2 * (1 - pnorm(boot.res$zs))
  return(list(boot.res, holdbeta))  
}

# Block Bootstrap with Controls
blockbootc <- function(df, sims = 1000, model, clustername, status, wt = "noweight"){
  holdbeta <- NULL
  
  for(s in 1:sims){
    # Choose clusters to sample
    
    # Make sure we get correct number of preclearance and non-preclearance states -
    # there are going to be some states that have preclearance counties within them but these might be sampled in the uncovered states, with the covered counties still marked as 1
    
    comp <- data.frame(table(df[, clustername], status))
    
    # Identify States in Treatment and Control Groups; Ultimately We Want Balance Here  
    t.clusters <- as.character(comp$Var1[comp$status == 1 & comp$Freq !=0])
    c.clusters <- as.character(unique(comp$Var1[!(comp$Var1 %in% t.clusters)]))
    
    # Get a Sample with Repeats  
    samp.clusters <- data.frame(state_name = as.character(sort(c(sample(t.clusters, length(t.clusters), replace = T), sample(c.clusters, length(c.clusters), replace = T)))))
    
    # Grab this sample out of the original data with repeats  
    bootdat <- apply(samp.clusters, 1, function(x) {df[which(df[ , clustername] == x),]})
    
    # Give every state a unique state ID. The reason to do this is that we want unique clusters. If I draw 3 Texases, for instance, I want 3 separate clusters since these are meant to be 
    # 3 distinct states despite having the same data
    for(i in 1:length(bootdat)){
      bootdat[[i]]$unitid <- paste(i, bootdat[[i]]$fips, sep = ".")
    }
    # Note that this generates a "unitid" field in the function. The model form is going to call a unitid variable from within the bootdat, so you'll want to have something
    #by this name in the bootstrapped data. This, as well as the samp.clusters line, also requires your df to have a column called "State_Name". We can always make the function more general to
    # let you specify what the cluster is, but I'll add that later if we need it.
    
    # Put the list into a single data frame  
    bootdat <- bootdat %>% bind_rows()
    
    # Run Regression on Data Using FELM  
    mod.form <- sub(".*= *(.*?) *=.*", "\\1", model$call)[2]
    if(wt == "noweight"){
      bootholder <- felm(as.formula(mod.form), data = bootdat)
    } else {
      bootholder <- felm(as.formula(mod.form), weights = bootdat[, wt], data = bootdat)
    }
    
    holdbeta <- rbind(holdbeta, bootholder$coefficients)
    print(s)
  }
  # Report Results
  boot.res <- data.frame(vars = unique(rownames(holdbeta)))
  boot.res$ses <- NA
  boot.res$orig.coef <- NA
  boot.res$orig.se <- NA
  for(p in 1:nrow(boot.res)){
    boot.res$ses[p] <- apply(as.matrix(holdbeta[row.names(holdbeta) == boot.res$vars[p], ]), 2, sd, na.rm = T)
    boot.res$orig.coef[p] <- coef(model)[names(coef(model)) == boot.res$vars[p]]
    boot.res$orig.se[p] <- summary(model)$coefficients[row.names(summary(model)$coefficients) == boot.res$vars[p], 2]
  }
  boot.res$lbounds <- boot.res$orig.coef - qnorm(0.975) * boot.res$se
  boot.res$ubounds <- boot.res$orig.coef + qnorm(0.975) * boot.res$se
  boot.res$zs <- abs(boot.res$orig.coef / boot.res$se)
  boot.res$ps <- 2 * (1 - pnorm(boot.res$zs))
  return(list(boot.res, holdbeta))  
}

# String Processing
splitit <- function(x,splitchar,n) sapply(strsplit(as.character(x), splitchar), "[", n)

trimit <- function(x) gsub("(^ +)|( +$)", "",x)

substrRight <- function(x, n){
  substr(x, nchar(x)-n+1, nchar(x))
}

# Factor to Numeric
as.numeric.factor <- function(x) {as.numeric(levels(x))[x]}

# Mode
getmode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}